Code
knitr::opts_chunk$set(
message = FALSE
)knitr::opts_chunk$set(
message = FALSE
)#if (!require("devtools", quietly = TRUE)) {
# install.packages("devtools")
#}
#devtools::install_github("adrientaudiere/MiscMetabar")library(MiscMetabar)
library(targets)
library(metacoder)
library(DT)
library(knitr)
library(ggplot2)
library(patchwork)
library(vegan)
library(FactoMineR)
library(factoextra)
library(DESeq2)
library(indicspecies)
library(sf)
library(adespatial)
library(ade4)
library(adegraphics)
library(kableExtra)
library(formattable)
library(microbiomeMarker)
library(ComplexUpset)d <- tar_read(d_vs_blast)d@sam_data$rank <- gsub("Weeding", "Grassing", d@sam_data$rank)
d@sam_data$inter_rank <- gsub("Weeding", "Grassing", d@sam_data$inter_rank)
d@sam_data$soil_practice <- gsub("Weeding", "Grassing", d@sam_data$soil_practice)
d@sam_data$terroir <- gsub("Perrier", "Vergèze", d@sam_data$terroir)d@sam_data$nb_spores <- rowMeans(cbind(as.numeric(d@sam_data$spores_rep1), as.numeric(d@sam_data$spores_rep2),as.numeric(d@sam_data$spores_rep3)),na.rm=TRUE)d@sam_data$organic <- ifelse(!d@sam_data$practice %in%
c("Conventional", "Conversion"),
"Organic", "NonOrganic"
)d@sam_data$paired_name <-
apply(as.matrix(d@sam_data[, c(5:16)]), 1, paste,
collapse = "--"
)
d@sam_data$paired_name[d@sam_data$terroir == "Vergèze"] <-
gsub("R", "", d@sam_data$ref_mycea[d@sam_data$terroir == "Vergèze"])tib_sam_data <- as_tibble(d@sam_data) %>%
mutate(across(starts_with("Myc_freq"), as.numeric)) %>%
mutate(across(starts_with("Myc_intensity_colonization"), as.numeric)) %>%
mutate(across(starts_with("Myc_intensity_root"), as.numeric)) %>%
mutate(across(starts_with("Arbuscul_richness"), as.numeric)) %>%
mutate(across(starts_with("Arbuscul_abundance"), as.numeric)) %>%
mutate(across(starts_with("Vesicle_richness"), as.numeric)) %>%
mutate(across(starts_with("Vesicle_abundance"), as.numeric)) %>%
mutate(Myc_freq = rowMeans(select(., starts_with("Myc_freq")), na.rm = TRUE)) %>%
mutate(Myc_intensity_colonization = rowMeans(select(
., starts_with("Myc_intensity_colonization")
), na.rm = TRUE)) %>%
mutate(Myc_intensity_root = rowMeans(select(., starts_with(
"Myc_intensity_root"
)), na.rm = TRUE)) %>%
mutate(Arbuscul_richness = rowMeans(select(., starts_with(
"Arbuscul_richness"
)), na.rm = TRUE)) %>%
mutate(Arbuscul_abundance = rowMeans(select(., starts_with(
"Arbuscul_abundance"
)), na.rm = TRUE)) %>%
mutate(Vesicle_richness = rowMeans(
select(., starts_with(
"Vesicle_richness"
)),
na.rm = TRUE
)) %>%
mutate(Vesicle_abundance = rowMeans(
select(., starts_with(
"Vesicle_abundance"
)),
na.rm = TRUE
)) %>% tibble::column_to_rownames(var = "X")
sample_data(d) <- tib_sam_datad <- clean_pq(subset_samples(d, compartment != "Mycelium"))
d_vs <- clean_pq(subset_samples(tar_read(d_vs), compartment != "Mycelium"))summary_plot_pq(d) +
labs(title = "Summary of final dataset") +
theme(plot.title = element_text(hjust = 0.5, size = 20, color = "#1e2b4c"))kable(tar_read(track_sequences_samples_clusters))| nb_sequences | nb_clusters | nb_samples | |
|---|---|---|---|
| Paired sequences | 8283716 | 28637 | 194 |
| Paired sequences without chimera | 7777304 | 9791 | 194 |
| Paired sequences without chimera and longer than 200bp | 7736170 | 7522 | 194 |
| ASV denoising | 7533009 | 7522 | 186 |
| OTU after vsearch reclustering at 97% | 7533009 | 1081 | 186 |
| OTU after blast filter without reclustering | 3819912 | 3793 | 182 |
| OTU after blast filter with reclustering | 3822304 | 219 | 182 |
d_muco <- clean_pq(subset_taxa(d, Class_PR2 == "Mucoromycota"))
d_vs_blast <- clean_pq(subset_taxa(tar_read(d_vs_blast), Class_PR2 == "Mucoromycota"))for(i in c(19:46)){
d_muco@sam_data[,i] <- as.numeric(unlist(d_muco@sam_data[,i]))
}d_roots <- clean_pq(subset_samples(d_muco, compartment == "Roots"))
d_spores <- clean_pq(subset_samples(d_muco, compartment == "Spores"))d_rs_merged_samples <- clean_pq(merge_samples2(
d_muco,
"paired_name",
default_fun =
function(x){MiscMetabar::diff_fct_diff_class(x, na.rm = T)}
))d_rs <- clean_pq(subset_samples_pq(
d_rs_merged_samples,
sample_sums(d_rs_merged_samples) > 2000
))# Figure S8
multitax_bar_pq(d_vs, "Supergroup_PR2",
"Subdivision_PR2", "Class_PR2",
nb_seq = TRUE
) +
ggtitle("Number of sequences (log10) including non-Mucoromycota")# Figure S6
dir.create("output/krona")Warning in dir.create("output/krona"): 'output/krona' existe déjà
krona(
d_vs,
ranks = c(11:19),
"output/krona/OTU_post_clustered_allEuk",
name = "OTU_post_clustered_allEuk"
)
krona(
d,
ranks = c(11:19),
"output/krona/OTU_filtOnMaarjam_AM",
name = "OTU_filtOnMaarjam_AM"
)
krona(
d_muco,
ranks = c(11:19),
"output/krona/OTU_filtOnMaarjamAndPr2_AM",
name = "OTU_filtOnMaarjamAndPr2_AM"
)
merge_krona(
c(
"output/krona/OTU_post_clustered_allEuk",
"output/krona/OTU_filtOnMaarjam_AM",
"output/krona/OTU_filtOnMaarjamAndPr2_AM"
),
"output/krona/Euk_to_AM_filtering_nb_seq.html"
)# Figure S7
krona(
d_vs,
ranks = c(11:19),
"output/krona/OTU_post_clustered_allEuk_Nbotu",
name = "OTU_post_clustered_allEuk",
nb_seq = F
)
krona(
d,
ranks = c(11:19),
"output/krona/OTU_filtOnMaarjam_AM_Nbotu",
name = "OTU_filtOnMaarjam_AM",
nb_seq = F
)
krona(
d_muco,
ranks = c(11:19),
"output/krona/OTU_filtOnMaarjamAndPr2_AM_Nbotu",
name = "OTU_filtOnMaarjamAndPr2_AM",
nb_seq = F
)
merge_krona(
c(
"output/krona/OTU_post_clustered_allEuk_Nbotu",
"output/krona/OTU_filtOnMaarjam_AM_Nbotu",
"output/krona/OTU_filtOnMaarjamAndPr2_AM_Nbotu"
),
"output/krona/Euk_to_AM_filtering_Nbotu.html"
)psm_otu <- psmelt(as_binary_otu_table(d_spores)) |>
filter(Abundance > 0) |>
group_by(Sample) |>
summarise(
"Abundance" = sum(Abundance),
"region" = region[1],
"nb_spores" = nb_spores[1]
)psm_res <- psmelt_samples_pq(d_spores) %>%
mutate(nb_spores_log10 = log10(nb_spores))# Fig S2
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_spores, rngseed = 221013)) %>% mutate(nb_spores_log10 = log10(nb_spores))
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Abundance, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_0, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_1, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_2, type = "non-parametric")# Figure S5
ggstatsplot::ggbetweenstats(psm_res, practice, nb_spores_log10, type = "non-parametric") # Figure S4
ggstatsplot::ggbetweenstats(
mutate(psm_res,
terroir = reorder(psm_res$terroir, psm_res$nb_spores)
),
terroir,
nb_spores_log10,
type = "non-parametric",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist"
) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Myc_freq, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Myc_intensity_colonization, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Arbuscul_abundance, type = "non-parametric")psm_res <- psmelt_samples_pq(d_rs)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_rs, rngseed = 221013))
# Figure S3a
(ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_2, type = "non-parametrique")) /
(ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_2, type = "non-parametrique"))# Figure S3b
(ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_2, type = "non-parametrique")) /
(ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_2, type = "non-parametrique"))# Figure S3c
(ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_2, type = "non-parametrique")) /
(ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_2, type = "non-parametrique"))ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Myc_freq, type = "non-parametrique", centrality.plotting = F) +
ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Myc_intensity_colonization, type = "non-parametrique", centrality.plotting = F) +
ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Arbuscul_abundance, type = "non-parametrique", centrality.plotting = F)ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Myc_freq, type = "non-parametrique", centrality.plotting = F) +
ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Myc_intensity_colonization, type = "non-parametrique", centrality.plotting = F) +
ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Arbuscul_abundance, type = "non-parametrique", centrality.plotting = F)psm_res <- psmelt_samples_pq(d_muco)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_muco, rngseed = 221013))ggstatsplot::ggbetweenstats(psm_res_rarefied,
compartment,
Hill_0,
type = "non-parametrique"
) +
ggstatsplot::ggbetweenstats(psm_res_rarefied,
compartment,
Hill_1,
type = "non-parametrique"
) +
ggstatsplot::ggbetweenstats(psm_res_rarefied,
compartment,
Hill_2,
type = "non-parametrique"
)ggplot(psm_res_rarefied, aes(y=Hill_0, x=compartment, color=practice))+
geom_point() +
geom_line(aes(group = paired_name)) +
facet_wrap(~terroir)psm_res <- psmelt_samples_pq(d_rs)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_rs, rngseed = 221013))# Figure 4
(ggstatsplot::ggbetweenstats(
mutate(
psm_res_rarefied,
terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_0)
),
terroir,
Hill_0,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_0))),levels(as.factor(psm_res_rarefied$terroir)))]) + coord_flip()) /
(ggstatsplot::ggbetweenstats(
mutate(
psm_res_rarefied,
terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_1)
),
terroir,
Hill_1,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_1))),levels(as.factor(psm_res_rarefied$terroir)))]) + coord_flip()) /
(ggstatsplot::ggbetweenstats(
mutate(
psm_res_rarefied,
terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_2)
),
terroir,
Hill_2,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_2))),levels(as.factor(psm_res_rarefied$terroir)))]) + coord_flip()) + plot_layout(axes = "collect")# Figure S10
(ggstatsplot::ggbetweenstats(
mutate(
psm_res,
terroir = reorder(psm_res$terroir, psm_res$Hill_0)
),
terroir,
Hill_0,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_0))),levels(as.factor(psm_res$terroir)))]) + coord_flip()) /
(ggstatsplot::ggbetweenstats(
mutate(
psm_res,
terroir = reorder(psm_res$terroir, psm_res$Hill_1)
),
terroir,
Hill_1,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_1))),levels(as.factor(psm_res$terroir)))]) + coord_flip()) /
(ggstatsplot::ggbetweenstats(
mutate(
psm_res,
terroir = reorder(psm_res$terroir, psm_res$Hill_2)
),
terroir,
Hill_2,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_2))),levels(as.factor(psm_res$terroir)))]) + coord_flip()) + plot_layout(axes = "collect")# Fig_7a
p_accu_balanced_mod_terroir <-
MiscMetabar::accu_plot_balanced_modality(d_rs,
"terroir",
nperm = 999,
step = 300,
quantile_prob=0.9,
sample.size = 10000,
progress_bar = FALSE)Warning in vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE): most
observed count data have counts 1, but smallest count is 3
Warning in vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE): most
observed count data have counts 1, but smallest count is 3
Warning in vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE): most
observed count data have counts 1, but smallest count is 3
Warning in vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE): most
observed count data have counts 1, but smallest count is 17
Warning in max(dff$xlab, na.rm = TRUE): aucun argument pour max ; -Inf est
renvoyé
p_accu_balanced_mod_terroir +
xlab("Number of sequences") +
ylab("Number of OTUs") +
ggrepel::geom_label_repel(data=summarise(group_by(p_accu_balanced_mod_terroir$data,factor) ,y=max(X1,na.rm = T)),aes(label=factor,y=y,x=max(p_accu_balanced_mod_terroir$data$x,na.rm = T)), max.overlaps = 1900, force = 20) + xlim(c(1,25000)) + theme(legend.position = "none" )Warning: Removed 14 rows containing missing values or values outside the scale range
(`geom_line()`).
# Figure 5
ggstatsplot::ggbetweenstats(psm_res_rarefied,
practice,
Hill_0,
centrality.plotting = F,
type = "non-parametrique"
) + ylab("Hill 0 (richness)") +
ggstatsplot::ggbetweenstats(psm_res_rarefied,
practice,
Hill_1,
centrality.plotting = F,
type = "non-parametrique"
) + ylab("Hill 1 (~Shannon)") +
ggstatsplot::ggbetweenstats(psm_res_rarefied,
practice,
Hill_2,
centrality.plotting = F,
type = "non-parametrique"
) + ylab("Hill 2 (~Simpson)") + plot_layout(axes = "collect")# Figure S11
ggstatsplot::ggbetweenstats(psm_res,
practice,
Hill_0,
centrality.plotting = F,
type = "non-parametrique"
) +
ggstatsplot::ggbetweenstats(psm_res,
practice,
Hill_1,
centrality.plotting = F,
type = "non-parametrique"
) +
ggstatsplot::ggbetweenstats(psm_res,
practice,
Hill_2,
centrality.plotting = F,
type = "non-parametrique"
)# Fig_7b
p_accu_balanced_mod <-
MiscMetabar::accu_plot_balanced_modality(d_rs,
"practice",
nperm = 999,
step = 300,
quantile_prob=0.90,
progress_bar = FALSE)Warning in max(dff$xlab, na.rm = TRUE): aucun argument pour max ; -Inf est
renvoyé
p_accu_balanced_mod +
xlab("Number of sequences") +
ylab("Number of OTUs") +
ggrepel::geom_label_repel(data=summarise(group_by(p_accu_balanced_mod$data,factor) ,y=max(X1,na.rm = T)),aes(label=factor,y=y,x=max(p_accu_balanced_mod$data$x,na.rm = T)), max.overlaps = 1900, force = 20) + xlim(c(1,18000)) + theme(legend.position = "none" )Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_line()`).
otu_hill <- vegan::renyi(d_rs@otu_table, scale = c(0, 1, 2), hill = TRUE)
low_div_samples <-
rownames(otu_hill[otu_hill$`0` < 10 |
otu_hill$`1` < 2 |
otu_hill$`2` < 2, ])
low_div_samples [1] "2019--Languedoc--Cevennes--St Maurice de Cazevieille--Conventional--Conventional--4.247388--44.035228--Conventional--Herbicide--Scratching--Herbicide"
[2] "PE13"
[3] "2023--Champagne--Côte des Blancs--Chetillons--Chetillons Desh sous rang_plein chimique--Conventional--4.034283--48.945738--Conventional--Herbicide--Herbicide--Herbicide"
[4] "2023--Champagne--Côte des Blancs--Chetillons--Chetillons Ecorces--Conventional--4.034333--48.945724--Conventional--Herbicide--Grassing--Herbicide"
[5] "2023--Champagne--Côte des Blancs--Cramant--BIONNES JEB--Conversion--4.018985--48.988723--Conversion--Scratching--Grassing--Grassing"
[6] "2023--Champagne--Côte des Blancs--Cramant--BATEAU AN39 JEB--Conversion--4.021354--48.995715--Conversion--Scratching--Grassing--Grassing"
[7] "2023--Champagne--Côte des Blancs--Cramant--TDB AM29 JEB--Conversion--4.016860--48.996831--Conversion--Scratching--Grassing--Grassing"
[8] "2023--Champagne--Côte des Blancs--Mireaux--JEB VOISIN 1 MIREAUX DESSOUS--Conventional--3.994335--48.984188--Conventional--Herbicide--Herbicide--Herbicide"
[9] "2023--Champagne--Côte des Blancs--Avize--AVIZE Desh chimique en plein--Conventional--4.014914--48.976412--Conventional--Herbicide--Herbicide--Herbicide"
[10] "2023--Champagne--Côte des Blancs--Avize--AVIZE Travail meca sous rang--Organic--4.014914--48.976351--Organic grassed--Scratching--Grassing--Grassing"
[11] "2020--Champagne--Côte des Bar--Les Riceys--Conventional--Conventional--4.392132--48.009899--Conventional--Herbicide--Scratching--Herbicide"
[12] "2020--Bordelais--Langoiran--Langoiran--Conventional chimique--Conventional---0.356503--44.713204--Conventional--Herbicide--Herbicide--Herbicide"
[13] "2020--Champagne--Vallee de la Marne--Vincelles--Conventional--Conventional--3.931795--49.077467--Conventional--Herbicide--Scratching--Herbicide"
[14] "2020--Camargue--Montcalm--Montcalm--Conventional--Conventional--4.320967--43.584075--Conventional--Scratching--Scratching--Scratching"
sum(otu_hill$`0` < 10)[1] 11
sum(otu_hill$`1` < 2)[1] 8
sum(otu_hill$`2` < 2)[1] 11
length(low_div_samples)[1] 14
low_div_samples_table <- tibble(cbind(as_tibble(d_rs@sam_data[low_div_samples, c("terroir", "practice", "rank", "inter_rank", "compartment")]), apply(otu_hill[low_div_samples, ], 2, round, digits = 2), sample_sums(d_rs)[low_div_samples]))Warning in class(x) <- tibble_class: La class(x) est une chaîne multiple
("tbl_df", "tbl", ...) ; le résultat ne sera plus un objet S4
low_div_samples_table$compartment <- ifelse(is.na(low_div_samples_table$compartment), "Spores + Roots", "Spores")
colnames(low_div_samples_table) <- c("Terroir", "Global practice", "Rank practice", "Inter rank practice", "Compartment", "Hill 0 (richness)", "Hill 1 (~Shannon)", "Hill 2 (~Simpson)", "nb_seq")
low_div_samples_table <- arrange(low_div_samples_table, as.numeric(`Hill 0 (richness)`))
low_div_samples_table <-
rbind(
low_div_samples_table,
c(
"",
"",
"",
"MEAN",
"(low -diversity samples)",
round(mean(low_div_samples_table$`Hill 0 (richness)`), 2),
round(mean(low_div_samples_table$`Hill 1 (~Shannon)`), 2),
round(mean(low_div_samples_table$`Hill 2 (~Simpson)`), 2),
round(mean(low_div_samples_table$nb_seq), 2)
),
c(
"",
"",
"",
"MEAN",
"(all samples)",
round(mean(otu_hill$`0`), 2),
round(mean(otu_hill$`1`), 2),
round(mean(otu_hill$`2`), 2),
round(mean(sample_sums(d_rs)), 2)
),
c(
"",
"",
"",
"MAX",
"(all samples)",
round(max(otu_hill$`0`), 2),
round(max(otu_hill$`1`), 2),
round(max(otu_hill$`2`), 2),
round(max(sample_sums(d_rs)), 2)
)
)
# Table 3
kbl(low_div_samples_table) |>
kable_classic(full_width = F, html_font = "Cambria")| Terroir | Global practice | Rank practice | Inter rank practice | Compartment | Hill 0 (richness) | Hill 1 (~Shannon) | Hill 2 (~Simpson) | nb_seq |
|---|---|---|---|---|---|---|---|---|
| Côte des Blancs | Conversion | Scratching | Grassing | Spores + Roots | 3 | 1.87 | 1.68 | 42775 |
| Côte des Blancs | Conventional | Herbicide | Grassing | Spores + Roots | 4 | 1.71 | 1.39 | 28481 |
| Langoiran | Conventional | Herbicide | Herbicide | Spores + Roots | 4 | 2.38 | 1.87 | 12643 |
| Côte des Blancs | Organic | Scratching | Grassing | Spores + Roots | 5 | 1.99 | 1.49 | 35728 |
| Vergèze | Organic | Scratching | Ploughing | Spores | 6 | 4.33 | 3.65 | 6768 |
| Côte des Blancs | Conventional | Herbicide | Herbicide | Spores + Roots | 6 | 1.34 | 1.12 | 41991 |
| Côte des Blancs | Conversion | Scratching | Grassing | Spores + Roots | 7 | 2.74 | 2.13 | 49313 |
| Côte des Blancs | Conventional | Herbicide | Herbicide | Spores + Roots | 8 | 1.68 | 1.26 | 46114 |
| Côte des Blancs | Conversion | Scratching | Grassing | Spores + Roots | 8 | 1.17 | 1.05 | 45845 |
| Vallee de la Marne | Conventional | Herbicide | Scratching | Spores + Roots | 9 | 2.7 | 1.78 | 11777 |
| Montcalm | Conventional | Scratching | Scratching | Spores | 9 | 3.45 | 2.94 | 24675 |
| Côte des Blancs | Conventional | Herbicide | Herbicide | Spores + Roots | 10 | 1.45 | 1.16 | 45017 |
| Cevennes | Conventional | Herbicide | Scratching | Spores + Roots | 12 | 1.91 | 1.35 | 93267 |
| Côte des Bar | Conventional | Herbicide | Scratching | Spores + Roots | 17 | 2.43 | 1.51 | 28971 |
| MEAN | (low -diversity samples) | 7.71 | 2.22 | 1.74 | 36668.93 | |||
| MEAN | (all samples) | 19.87 | 6.23 | 4.35 | 48209.53 | |||
| MAX | (all samples) | 47 | 15.06 | 9.04 | 101104 |
# Figure S9
p <- phyloseq::plot_ordination(d_muco,
vegan::decorana(vegdist(as(otu_table(d_muco), "matrix"),
method = "robust.aitchison"
)),
color = "region",
shape = "compartment"
) +
geom_point(size = 3) +
stat_ellipse(inherit.aes = F, aes(x = DCA1, y = DCA2, linetype = compartment))Warning in phyloseq::plot_ordination(d_muco,
vegan::decorana(vegdist(as(otu_table(d_muco), : `Ordination species/OTU/taxa
coordinate indices did not match `physeq` index names. Setting corresponding
coordinates to NULL.
p + geom_line(data = p$data, aes(group = paired_name))p <-
phyloseq::plot_ordination(d_rs,
ordinate(d_rs,
method = "NMDS",
distance = "bray"
),
color = "terroir",
shape = "practice"
) + facet_wrap("region") + scale_shape_manual(values=c(25,23,24))Square root transformation
Wisconsin double standardization
Run 0 stress 0.2485558
Run 1 stress 0.2596844
Run 2 stress 0.2663072
Run 3 stress 0.2567211
Run 4 stress 0.256791
Run 5 stress 0.2626177
Run 6 stress 0.2636841
Run 7 stress 0.2612171
Run 8 stress 0.2652946
Run 9 stress 0.262289
Run 10 stress 0.2523773
Run 11 stress 0.2583239
Run 12 stress 0.2599351
Run 13 stress 0.2562442
Run 14 stress 0.409148
Run 15 stress 0.2593146
Run 16 stress 0.2658447
Run 17 stress 0.2580753
Run 18 stress 0.2626311
Run 19 stress 0.261243
Run 20 stress 0.2584539
*** Best solution was not repeated -- monoMDS stopping criteria:
2: no. of iterations >= maxit
18: stress ratio > sratmax
p_nmds <- p + geom_point(
data = select(p$data, -c(region)),
fill = "grey40",
color = "grey20",
size = 1,
alpha = 0.5
) +
geom_point(size = 3, aes(fill=terroir)) +
ggtitle("NMDS on bray distance") +
scale_fill_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0")) +
scale_color_manual(values=rep("grey20",14))
# Figure 9
p_nmds# Figure 6a
upset_pq(d_rs, "practice", taxa_fill = "Family",
set_sizes=(
upset_set_size()
+ geom_text(aes(label=after_stat(count)), hjust=-0.3, color="white", stat='count')
))ggvenn_pq(d_rs, "practice")# Figure 6b
upset_pq(d_rs, "terroir", taxa_fill = "Family", min_size = 2, height_ratio = 0.6
,set_sizes=(
upset_set_size()
+ geom_label(aes(label=after_stat(count)), hjust=-0.3, size=2.5, stat='count')
))Warning: Removed 61 rows containing non-finite outside the scale range
(`stat_count()`).
library(geodist)
lat_lon <- as_tibble(d_rs@sam_data[, c("lat", "long")])
lat_lon$lat <- as.numeric(lat_lon$lat)
lat_lon$long <- as.numeric(lat_lon$long)
dist_spatial_meter <- as.dist(geodist(lat_lon, measure = "geodesic"),
upper = FALSE
)lon_lat_rs <- d_rs@sam_data[, c("long", "lat")]
lon_lat_rs$long <- as.numeric(lon_lat_rs$long)
lon_lat_rs$lat <- as.numeric(lon_lat_rs$lat)
MEM <- dbmem(lon_lat_rs, MEM.autocor = "non-null")test_MEM <- moran.randtest(MEM, nrepet = 999)
test_MEM$pvalue_adjust <- p.adjust(test_MEM$pvalue, method = "BH")d_rs@sam_data$MEM_1 <- MEM[, 1]
d_rs@sam_data$MEM_2 <- MEM[, 2]
res_ado_spatial_robAit <- adonis_pq(
d_rs,
"MEM_1 + MEM_2 + practice + inter_rank + rank + terroir",
correction_for_sample_size = TRUE,
dist_method = "robust.aitchison"
)
res_ado_spatial_robAit_rarefy <- adonis_pq(
rarefy_even_depth(d_rs, rngseed = 626),
"MEM_1 + MEM_2 + practice + inter_rank + rank + terroir",
correction_for_sample_size = FALSE,
dist_method = "robust.aitchison"
)
res_ado_spatial_bray <- adonis_pq(
d_rs,
"MEM_1 + MEM_2 + practice + inter_rank + rank + terroir",
correction_for_sample_size = TRUE,
dist_method = "bray"
)
res_ado_spatial_bray_rarefy <- adonis_pq(
rarefy_even_depth(d_rs, rngseed = 626),
"MEM_1 + MEM_2 + practice + inter_rank + rank + terroir",
correction_for_sample_size = FALSE,
dist_method = "bray"
)
df <- data.frame(res_ado_spatial_bray)
df$names <-
factor(
rownames(df),
levels = c(
"sample_size",
"MEM_1",
"MEM_2",
"practice",
"inter_rank",
"rank",
"terroir",
"Residual",
"Total"
)
)
df <- df %>% filter(!names == "Total")
p1 <-
ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c() +
geom_text(
label = case_when(
df$Pr..F. < 0.001 ~ "***",
df$Pr..F. < 0.01 ~ "**",
df$Pr..F. < 0.05 ~ "*",
df$Pr..F. > 0.05 ~ "ns"
),
nudge_y = 0.01
) +
ggtitle("Permanova on Bray")
df <- data.frame(res_ado_spatial_robAit)
df$names <-
factor(
rownames(df),
levels = c(
"MEM_1",
"MEM_2",
"practice",
"inter_rank",
"rank",
"terroir",
"Residual",
"Total"
)
)
df <- df %>% filter(!names == "Total")
p2 <-
ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c() +
geom_text(
label = case_when(
df$Pr..F. < 0.001 ~ "***",
df$Pr..F. < 0.01 ~ "**",
df$Pr..F. < 0.05 ~ "*",
df$Pr..F. > 0.05 ~ "ns"
),
nudge_y = 0.01
) +
ggtitle("Permanova on robust Aitchison")
df <- data.frame(res_ado_spatial_bray_rarefy)
df$names <-
factor(
rownames(df),
levels = c(
"MEM_1",
"MEM_2",
"practice",
"inter_rank",
"rank",
"terroir",
"Residual",
"Total"
)
)
df <- df %>% filter(!names == "Total")
p3 <-
ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c() +
geom_text(
label = case_when(
df$Pr..F. < 0.001 ~ "***",
df$Pr..F. < 0.01 ~ "**",
df$Pr..F. < 0.05 ~ "*",
df$Pr..F. > 0.05 ~ "ns"
),
nudge_y = 0.01
) +
ggtitle("Permanova on Bray after rarefaction")
df <- data.frame(res_ado_spatial_robAit_rarefy)
df$names <-
factor(
rownames(df),
levels = c(
"MEM_1",
"MEM_2",
"practice",
"inter_rank",
"rank",
"terroir",
"Residual",
"Total"
)
)
df <- df %>% filter(!names == "Total")
p4 <-
ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c() +
geom_text(
label = case_when(
df$Pr..F. < 0.001 ~ "***",
df$Pr..F. < 0.01 ~ "**",
df$Pr..F. < 0.05 ~ "*",
df$Pr..F. > 0.05 ~ "ns"
),
nudge_y = 0.01
) +
ggtitle("Permanova on robust Aitchison after rarefaction")
(p1 + ylim(0, 0.12) + p2 + ylim(0, 0.12)) /
(p3 + ylim(0, 0.12) + p4 + ylim(0, 0.12))Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), d_rs@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 0.16217 0.012475 0.6941 0.7615
Residuals 61 1.09633 0.017973
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), d_rs@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.05732 0.028658 2.1729 0.1213
Residuals 72 0.94962 0.013189
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 0.16217 0.012475 0.6941 0.7615
Residuals 61 1.09633 0.017973
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.05732 0.028658 2.1729 0.1213
Residuals 72 0.94962 0.013189
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 393.8 30.2920 9.5001 2.306e-10 ***
Residuals 61 194.5 3.1886
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 39.37 19.685 1.4366 0.2445
Residuals 72 986.60 13.703
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 393.8 30.2920 9.5001 2.306e-10 ***
Residuals 61 194.5 3.1886
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 39.37 19.685 1.4366 0.2445
Residuals 72 986.60 13.703
soil_prop <- as_tibble(d_rs@sam_data) |>
select(paired_name, practice, organic, terroir, Coarse_sand:CEC)
soil_prop_num <- soil_prop |>
select(-all_of(c("practice", "organic", "terroir", "CaO", "Total_sand", "Total_filt"))) |>
tibble::column_to_rownames("paired_name") |>
mutate(across(everything(), as.numeric)) |>
tidyr::drop_na()pca_test_res <- PCAtest::PCAtest(soil_prop_num, varcorr = T, plot = F, counter = F)
Sampling bootstrap replicates... Please wait
Calculating confidence intervals of empirical statistics... Please wait
Sampling random permutations... Please wait
Comparing empirical statistics with their null distributions... Please wait
========================================================
Test of PCA significance: 17 variables, 44 observations
1000 bootstrap replicates, 1000 random permutations
========================================================
Empirical Psi = 55.8085, Max null Psi = 9.0868, Min null Psi = 3.6595, p-value = 0
Empirical Phi = 0.4530, Max null Phi = 0.1828, Min null Phi = 0.1160, p-value = 0
Empirical eigenvalue #1 = 7.28487, Max null eigenvalue = 2.88796, p-value = 0
Empirical eigenvalue #2 = 3.15651, Max null eigenvalue = 2.39696, p-value = 0
Empirical eigenvalue #3 = 2.59627, Max null eigenvalue = 2.09788, p-value = 0
Empirical eigenvalue #4 = 1.16046, Max null eigenvalue = 1.87236, p-value = 1
Empirical eigenvalue #5 = 0.93253, Max null eigenvalue = 1.59515, p-value = 1
Empirical eigenvalue #6 = 0.74202, Max null eigenvalue = 1.48039, p-value = 1
Empirical eigenvalue #7 = 0.39204, Max null eigenvalue = 1.32508, p-value = 1
Empirical eigenvalue #8 = 0.24218, Max null eigenvalue = 1.18435, p-value = 1
Empirical eigenvalue #9 = 0.19075, Max null eigenvalue = 1.0485, p-value = 1
Empirical eigenvalue #10 = 0.08962, Max null eigenvalue = 0.99123, p-value = 1
Empirical eigenvalue #11 = 0.07595, Max null eigenvalue = 0.85026, p-value = 1
Empirical eigenvalue #12 = 0.05693, Max null eigenvalue = 0.7427, p-value = 1
Empirical eigenvalue #13 = 0.0369, Max null eigenvalue = 0.65767, p-value = 1
Empirical eigenvalue #14 = 0.02507, Max null eigenvalue = 0.55583, p-value = 1
Empirical eigenvalue #15 = 0.01767, Max null eigenvalue = 0.48237, p-value = 1
Empirical eigenvalue #16 = 0.00018, Max null eigenvalue = 0.40241, p-value = 1
Empirical eigenvalue #17 = 5e-05, Max null eigenvalue = 0.35079, p-value = 1
PC 1 is significant and accounts for 42.9% (95%-CI:38.8-49.7) of the total variation
PC 2 is significant and accounts for 18.6% (95%-CI:16.6-25) of the total variation
PC 3 is significant and accounts for 15.3% (95%-CI:10.5-17.9) of the total variation
The first 3 PC axes are significant and account for 76.7% of the total variation
Variables 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 15, 16, and 17 have significant loadings on PC 1
Variables 4, 8, 10, 12, and 14 have significant loadings on PC 2
Variables 5, 6, and 7 have significant loadings on PC 3
Variables 8, 10, 13, 15, and 17 have significant correlations with PC 1
Variables 8, 10, 12, and 14 have significant correlations with PC 2
Variables , and have significant correlations with PC 3
pval <- c()
for (i in seq_len(length(pca_test_res$`Percentage of variation of empirical PC's`))) {
obs <- pca_test_res$`Percentage of variation of empirical PC's`[i]
null_model <- pca_test_res$`Percentage of variation of randomized data`[, i]
pval[i] <- (sum(obs < null_model) + 1) / (1 + nrow(pca_test_res$`Percentage of variation of randomized data`))
}
pca_test_pval_adj <- p.adjust(pval, method = "BH")
pca_test_pval_adj [1] 0.005661006 0.005661006 0.005661006 1.000000000 1.000000000 1.000000000
[7] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
[13] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
res_pca <- dudi.pca(soil_prop_num, scannf = F, nf = 4)p_pca_variable <- fviz_pca_var(res_pca,
col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)
p_pca_variable_axe3_4 <- fviz_pca_var(res_pca,
col.var = "cos2", axes = c(3, 4),
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)
res_pca_var <- get_pca_var(res_pca)
# Fig. 2
(p_pca_variable + p_pca_variable_axe3_4) / free((factoextra::fviz_eig(res_pca)) +
(ggcorrplot::ggcorrplot(res_pca_var$cor)))practice <- unlist(lapply(strsplit(rownames(res_pca$tab), "--"), function(x) {
x[5]
}))
terroir <- unlist(lapply(strsplit(rownames(res_pca$tab), "--"), function(x) {
x[2]
}))
p_pca_ind_terroir <- fviz_pca_ind(res_pca,
habillage = terroir,
label = FALSE,
repel = TRUE,
mean.point = FALSE,
addEllipses = TRUE,
pointsize = 3,
ellipse.type = "convex"
) +
scale_shape_manual(values = rep(16, length(unique(terroir))))
p_pca_ind_terroir$data$terroir <- terroir
p_pca_ind_terroir$data$practice <- practice
p_pca_ind_terroir_practice <- p_pca_ind_terroir + geom_point(aes(x = x, y = y),
shape = case_when(
practice == "Organic" ~ 2,
practice == "Conventional" ~ 6,
practice == "Conversion" ~ 5
), size = 3.5,
color = rgb(0, 0, 0, 0.7)
)
# value
kw1_terroir <- kruskal.test(res_pca$tab[, 1], terroir)
kw2_terroir <- kruskal.test(res_pca$tab[, 2], terroir)
kw3_terroir <- kruskal.test(res_pca$tab[, 3], terroir)
kw1_practice <- kruskal.test(res_pca$tab[, 1], practice)
kw2_practice <- kruskal.test(res_pca$tab[, 2], practice)
kw3_practice <- kruskal.test(res_pca$tab[, 3], practice)
# Figure 3
p_pca_ind_terroir_practice +
annotate("text", x = 2, y = 4.5, label = paste("Kruskal-Wallis (axe1 - terroir): pval=", round(as.numeric(kw1_terroir$p.value),4)), size = 3) +
annotate("text", x = 2, y = 4.25, label = paste("Kruskal-Wallis (axe2 - terroir): pval=", round(as.numeric(kw2_terroir$p.value),4)), size = 3) +
annotate("text", x = 2, y = 4, label = paste("Kruskal-Wallis (axe3 - terroir): pval=", round(as.numeric(kw3_terroir$p.value),4)), size = 3) +
annotate("text", x = -2, y = 4.5, label = paste("Kruskal-Wallis (axe1 - practice): pval=", round(as.numeric(kw1_practice$p.value),4)), size = 3) +
annotate("text", x = -2, y = 4.25, label = paste("Kruskal-Wallis (axe2 - practice): pval=", round(as.numeric(kw2_practice$p.value),4)), size = 3) +
annotate("text", x = -2, y = 4, label = paste("Kruskal-Wallis (axe3 - practice): pval=", round(as.numeric(kw3_practice$p.value),4)), size = 3)Warning: Removed 21 rows containing missing values or values outside the scale range
(`geom_point()`).
# Figure S1
ggplot(d_rs@sam_data, aes(x = practice, y = as.numeric(Cu), fill = practice)) +
geom_violin() +
geom_jitter()Warning: Removed 23 rows containing non-finite outside the scale range
(`stat_ydensity()`).
Warning: Removed 23 rows containing missing values or values outside the scale range
(`geom_point()`).
kruskal.test(as.numeric(d_rs@sam_data$Cu), d_rs@sam_data$practice)
Kruskal-Wallis rank sum test
data: as.numeric(d_rs@sam_data$Cu) and d_rs@sam_data$practice
Kruskal-Wallis chi-squared = 0.72149, df = 2, p-value = 0.6972
d_with_pca <-
clean_pq(subset_samples_pq(d_rs, d_rs@sam_data$paired_name %in% rownames(soil_prop_num)))
res_pca_ind <- get_pca_ind(res_pca)
if (!identical(match(rownames(res_pca_ind$coord), d_with_pca@sam_data$paired_name), sort(match(rownames(res_pca_ind$coord), d_with_pca@sam_data$paired_name)))) {
stop("ERROR")
}
d_with_pca <- add_info_to_sam_data(d_with_pca, res_pca_ind$coord)dim_df <- as_tibble(d_with_pca@sam_data) |>
select(c(starts_with("Dim"), "nb_seq", "nb_otu"))
cor_results <- correlation::correlation(dim_df)
cor_results %>%
summary(redundant = TRUE) %>%
plot()lon_lat_with_pca <- d_with_pca@sam_data[, c("long", "lat")]
lon_lat_with_pca$long <- as.numeric(lon_lat_with_pca$long)
lon_lat_with_pca$lat <- as.numeric(lon_lat_with_pca$lat)
MEM_with_pca <- dbmem(lon_lat_with_pca, MEM.autocor = "non-null")
d_with_pca@sam_data$MEM_1 <- MEM_with_pca[, 1]
d_with_pca@sam_data$MEM_2 <- MEM_with_pca[, 2]res_ado_spatial_soil_bray <-
adonis_pq(
d_with_pca,
"MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
correction_for_sample_size = TRUE
)
res_ado_spatial_soil_robAit <-
adonis_pq(
d_with_pca,
"MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
correction_for_sample_size = TRUE,
dist_method = "robust.aitchison"
)
res_ado_spatial_soil_bray_rarefy <-
adonis_pq(
rarefy_even_depth(d_with_pca, rngseed = 626),
"MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
correction_for_sample_size = FALSE
)
res_ado_spatial_soil_robAit_rarefy <-
adonis_pq(
rarefy_even_depth(d_with_pca, rngseed = 626),
"MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
correction_for_sample_size = FALSE,
dist_method = "robust.aitchison"
)
anova(betadisper(phyloseq::distance(d_rs, method = "bray"), d_rs@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.05732 0.028658 2.1729 0.1213
Residuals 72 0.94962 0.013189
anova(betadisper(phyloseq::distance(d_rs, method = "bray"), d_rs@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 0.16217 0.012475 0.6941 0.7615
Residuals 61 1.09633 0.017973
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 39.37 19.685 1.4366 0.2445
Residuals 72 986.60 13.703
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 393.8 30.2920 9.5001 2.306e-10 ***
Residuals 61 194.5 3.1886
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# TABLE 2
kbl(res_ado_spatial_bray) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| sample_size | 1 | 1.0052204 | 0.0670420 | 7.319743 | 0.001 |
| MEM_1 | 1 | 1.5103278 | 0.1007296 | 10.997797 | 0.001 |
| MEM_2 | 1 | 0.9164972 | 0.0611247 | 6.673684 | 0.001 |
| practice | 2 | 0.4789872 | 0.0319455 | 1.743928 | 0.036 |
| inter_rank | 3 | 0.5442769 | 0.0362999 | 1.321092 | 0.128 |
| rank | 2 | 0.6046954 | 0.0403295 | 2.201614 | 0.004 |
| terroir | 13 | 2.9300449 | 0.1954160 | 1.641215 | 0.001 |
| Residual | 51 | 7.0038313 | 0.4671126 | NA | NA |
| Total | 74 | 14.9938812 | 1.0000000 | NA | NA |
kbl(res_ado_spatial_soil_bray) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| sample_size | 1 | 0.9095040 | 0.1022131 | 6.3231283 | 0.001 |
| MEM_1 | 1 | 0.2091359 | 0.0235034 | 1.4539720 | 0.149 |
| MEM_2 | 1 | 1.0896243 | 0.1224557 | 7.5753755 | 0.001 |
| Dim.1 | 1 | 0.3845781 | 0.0432202 | 2.6736956 | 0.007 |
| Dim.2 | 1 | 0.2413033 | 0.0271185 | 1.6776088 | 0.080 |
| Dim.3 | 1 | 0.0934249 | 0.0104994 | 0.6495160 | 0.783 |
| practice | 2 | 0.2815596 | 0.0316426 | 0.9787409 | 0.508 |
| inter_rank | 2 | 0.1453520 | 0.0163352 | 0.5052642 | 0.957 |
| rank | 2 | 0.4576684 | 0.0514343 | 1.5909201 | 0.070 |
| terroir | 11 | 2.2092074 | 0.2482782 | 1.3962758 | 0.022 |
| Residual | 20 | 2.8767531 | 0.3232993 | NA | NA |
| Total | 43 | 8.8981110 | 1.0000000 | NA | NA |
# TABLE S2
kbl(res_ado_spatial_bray_rarefy) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| MEM_1 | 1 | 1.5224735 | 0.1079197 | 11.273304 | 0.001 |
| MEM_2 | 1 | 1.1172045 | 0.0791924 | 8.272450 | 0.001 |
| practice | 2 | 0.5880731 | 0.0416852 | 2.177222 | 0.011 |
| inter_rank | 3 | 0.3890128 | 0.0275750 | 0.960161 | 0.478 |
| rank | 2 | 0.5861573 | 0.0415495 | 2.170130 | 0.016 |
| terroir | 13 | 2.8818779 | 0.2042804 | 1.641473 | 0.003 |
| Residual | 52 | 7.0226637 | 0.4977978 | NA | NA |
| Total | 74 | 14.1074628 | 1.0000000 | NA | NA |
# TABLE S3
kbl(res_ado_spatial_soil_bray_rarefy) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| MEM_1 | 1 | 0.3924837 | 0.0479272 | 2.8426837 | 0.014 |
| MEM_2 | 1 | 1.0983511 | 0.1341225 | 7.9551458 | 0.001 |
| Dim.1 | 1 | 0.4058745 | 0.0495624 | 2.9396712 | 0.008 |
| Dim.2 | 1 | 0.2866462 | 0.0350031 | 2.0761235 | 0.054 |
| Dim.3 | 1 | 0.1150841 | 0.0140532 | 0.8335323 | 0.581 |
| practice | 2 | 0.2932271 | 0.0358067 | 1.0618937 | 0.388 |
| inter_rank | 2 | 0.1439472 | 0.0175778 | 0.5212911 | 0.945 |
| rank | 2 | 0.4694158 | 0.0573216 | 1.6999441 | 0.061 |
| terroir | 11 | 2.0847081 | 0.2545690 | 1.3726490 | 0.056 |
| Residual | 21 | 2.8994281 | 0.3540566 | NA | NA |
| Total | 43 | 8.1891660 | 1.0000000 | NA | NA |
# TABLE S4
kbl(res_ado_spatial_robAit_rarefy) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| MEM_1 | 1 | 118.30117 | 0.0313656 | 2.7789327 | 0.001 |
| MEM_2 | 1 | 110.52944 | 0.0293050 | 2.5963722 | 0.002 |
| practice | 2 | 117.62397 | 0.0311860 | 1.3815125 | 0.043 |
| inter_rank | 3 | 93.00477 | 0.0246586 | 0.7282373 | 0.965 |
| rank | 2 | 154.47993 | 0.0409577 | 1.8143917 | 0.002 |
| terroir | 13 | 964.07334 | 0.2556078 | 1.7420303 | 0.001 |
| Residual | 52 | 2213.67750 | 0.5869192 | NA | NA |
| Total | 74 | 3771.69011 | 1.0000000 | NA | NA |
# TABLE S5
kbl(res_ado_spatial_soil_robAit_rarefy) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| MEM_1 | 1 | 112.72637 | 0.0339610 | 1.6387773 | 0.008 |
| MEM_2 | 1 | 112.45902 | 0.0338805 | 1.6348907 | 0.006 |
| Dim.1 | 1 | 191.81334 | 0.0577875 | 2.7885165 | 0.001 |
| Dim.2 | 1 | 93.30719 | 0.0281106 | 1.3564679 | 0.087 |
| Dim.3 | 1 | 112.44158 | 0.0338752 | 1.6346370 | 0.014 |
| practice | 2 | 128.09782 | 0.0385919 | 0.9311211 | 0.678 |
| inter_rank | 2 | 98.34228 | 0.0296275 | 0.7148331 | 0.956 |
| rank | 2 | 158.00090 | 0.0476008 | 1.1484814 | 0.175 |
| terroir | 11 | 867.57636 | 0.2613741 | 1.1465934 | 0.057 |
| Residual | 21 | 1444.52439 | 0.4351909 | NA | NA |
| Total | 43 | 3319.28925 | 1.0000000 | NA | NA |
# Figure 8a
res_varpart_rarefy <- var_par_rarperm_pq(
physeq = d_with_pca,
list_component = list(
"Spatial" = c("MEM_1", "MEM_2"),
"Soil" = c("Dim.1", "Dim.2", "Dim.3"),
"Practice" = c("practice", "inter_rank", "rank"),
"Terroir" = c("terroir")
),
nperm = 99,
dbrda_computation = TRUE,
progress_bar = FALSE
)
plot_var_part_pq(res_varpart_rarefy, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)# Figure 8b
res_varpart_rarefy_wo_soil <- var_par_rarperm_pq(
physeq = d_rs,
list_component = list(
"Spatial" = c("MEM_1", "MEM_2"),
"Practice" = c("practice", "inter_rank", "rank"),
"Terroir" = c("terroir")
),
nperm = 99,
dbrda_computation = TRUE,
progress_bar = FALSE
)
plot_var_part_pq(res_varpart_rarefy_wo_soil, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)res_varpart_robAit_rarefy <- var_par_rarperm_pq(
physeq = d_with_pca,
list_component = list(
"Spatial" = c("MEM_1", "MEM_2"),
"Soil" = c("Dim.1", "Dim.2", "Dim.3"),
"Practice" = c("practice", "inter_rank", "rank"),
"Terroir" = c("terroir")
),
dist_method = "robust.aitchison",
nperm = 99,
progress_bar = FALSE
)
plot_var_part_pq(res_varpart_robAit_rarefy, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)# Figure S12b
res_varpart_robAit_rarefy_wo_soil <- var_par_rarperm_pq(
physeq = d_rs,
list_component = list(
"Spatial" = c("MEM_1", "MEM_2"),
"Practice" = c("practice", "inter_rank", "rank"),
"Terroir" = c("terroir")
),
dist_method = "robust.aitchison",
nperm = 99,
dbrda_computation = TRUE,
progress_bar = FALSE
)
plot_var_part_pq(res_varpart_robAit_rarefy_wo_soil, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)# Figure 10
res_mpt_terroir <-
multipatt(as.matrix(d_rs@otu_table),
d_rs@sam_data$terroir,
control = how(nperm = 9999),
max.order = 3
)
res_mpt_terroir_df <- res_mpt_terroir$sign
res_mpt_terroir_df$p.adj <- p.adjust(res_mpt_terroir_df$p.value, method = "BH")
res_mpt_terroir_df$OTU_names <- rownames(res_mpt_terroir_df)
res_mpt_terroir_df_signif <-
res_mpt_terroir_df %>%
filter(p.adj < 0.05) %>%
tidyr::pivot_longer(cols = starts_with("s."))
tax_tab <- as.data.frame(d_rs@tax_table)
tax_tab$otu <- rownames(tax_tab)
res_mpt_terroir_df_signif_taxo <- left_join(res_mpt_terroir_df_signif,tax_tab, by = join_by("OTU_names" == "otu"))
ggplot(
res_mpt_terroir_df_signif_taxo,
aes(
x = OTU_names,
y = name,
alpha = value,
color = Genus
)
) +
geom_point(size = 6) +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
))+
scale_alpha(range = c(0, 1))res_mpt_terroir_rg <-
multipatt(as.matrix(rarefy_even_depth(d_rs, rngseed = 626)@otu_table),
d_rs@sam_data$terroir,
control = how(nperm = 9999),
max.order = 3,
func = "r.g"
)
res_mpt_terroir_rg_df <- res_mpt_terroir_rg$sign
res_mpt_terroir_rg_df$p.adj <- p.adjust(res_mpt_terroir_rg_df$p.value, method = "BH")
res_mpt_terroir_rg_df$OTU_names <- rownames(res_mpt_terroir_rg_df)
res_mpt_terroir_rg_df_signif <-
res_mpt_terroir_rg_df %>%
filter(p.adj < 0.05) %>%
tidyr::pivot_longer(cols = starts_with("s."))
ggplot(
res_mpt_terroir_rg_df_signif,
aes(
x = OTU_names,
y = name,
alpha = value,
color = stat
)
) +
geom_point(size = 4) +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
))res_mpt_practice <-
multipatt(as.matrix(d_rs@otu_table),
d_rs@sam_data$practice,
control = how(nperm = 9999),
max.order = 3
)
res_mpt_practice_df <- res_mpt_practice$sign
res_mpt_practice_df$p.adj <- p.adjust(res_mpt_practice_df$p.value, method = "BH")
res_mpt_practice_df$OTU_names <- rownames(res_mpt_practice_df)
res_mpt_practice_df_signif <-
res_mpt_practice_df %>%
filter(p.adj < 0.05) %>%
tidyr::pivot_longer(cols = starts_with("s."))
res_mpt_practice_df_signif# A tibble: 0 × 7
# ℹ 7 variables: index <int>, stat <dbl>, p.value <dbl>, p.adj <dbl>,
# OTU_names <chr>, name <chr>, value <dbl>
res_mpt_practice_rg <-
multipatt(as.matrix(d_rs@otu_table),
d_rs@sam_data$practice,
control = how(nperm = 9999),
max.order = 3,
func = "r.g"
)
res_mpt_practice_rg_df <- res_mpt_practice_rg$sign
res_mpt_practice_rg_df$p.adj <- p.adjust(res_mpt_practice_rg_df$p.value, method = "BH")
res_mpt_practice_rg_df$OTU_names <- rownames(res_mpt_practice_rg_df)
res_mpt_practice_rg_df_signif <-
res_mpt_practice_rg_df %>%
filter(p.adj < 0.05) %>%
tidyr::pivot_longer(cols = starts_with("s."))
res_mpt_practice_rg_df_signif# A tibble: 0 × 7
# ℹ 7 variables: index <int>, stat <dbl>, p.value <dbl>, p.adj <dbl>,
# OTU_names <chr>, name <chr>, value <dbl>
sessionInfo()R version 4.3.3 (2024-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 11 (bullseye)
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so; LAPACK version 3.9.0
locale:
[1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
[3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
[5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
[7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
time zone: Europe/Paris
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] geodist_0.0.8 ComplexUpset_1.3.3
[3] microbiomeMarker_1.8.0 formattable_0.2.1
[5] kableExtra_1.4.0 adegraphics_1.0-21
[7] ade4_1.7-22 adespatial_0.3-23
[9] sf_1.0-16 indicspecies_1.7.14
[11] DESeq2_1.42.1 SummarizedExperiment_1.32.0
[13] Biobase_2.62.0 MatrixGenerics_1.14.0
[15] matrixStats_1.3.0 GenomicRanges_1.54.1
[17] GenomeInfoDb_1.38.8 IRanges_2.36.0
[19] S4Vectors_0.40.2 BiocGenerics_0.48.1
[21] factoextra_1.0.7 FactoMineR_2.11
[23] vegan_2.6-6.1 lattice_0.22-6
[25] permute_0.9-7 patchwork_1.2.0
[27] knitr_1.46 DT_0.33
[29] metacoder_0.3.7 targets_1.7.0
[31] MiscMetabar_0.9.2 purrr_1.0.2
[33] dplyr_1.1.4 dada2_1.30.0
[35] Rcpp_1.0.12 ggplot2_3.5.1
[37] phyloseq_1.46.0
loaded via a namespace (and not attached):
[1] igraph_2.0.3 mia_1.10.0
[3] Formula_1.2-5 scater_1.30.1
[5] rematch2_2.1.2 zlibbioc_1.48.2
[7] tidyselect_1.2.1 bit_4.0.5
[9] doParallel_1.0.17 BWStest_0.2.3
[11] clue_0.3-65 rjson_0.2.21
[13] blob_1.2.4 stringr_1.5.1
[15] rngtools_1.5.2 S4Arrays_1.2.1
[17] parallel_4.3.3 RNeXML_2.4.11
[19] correlation_0.8.4 png_0.1-8
[21] cli_3.6.3 ggplotify_0.1.2
[23] CVXR_1.0-12 bayestestR_0.13.2
[25] multtest_2.58.0 MultiAssayExperiment_1.28.0
[27] bluster_1.12.0 BiocNeighbors_1.20.2
[29] TreeSummarizedExperiment_2.10.0 adegenet_2.1.10
[31] mime_0.12 evaluate_0.23
[33] tidytree_0.4.6 ComplexHeatmap_2.18.0
[35] ggh4x_0.2.8 stringi_1.8.4
[37] backports_1.4.1 PMCMRplus_1.9.10
[39] lmerTest_3.1-3 gsl_2.1-8
[41] XML_3.99-0.16.1 Exact_3.2
[43] ggVennDiagram_1.5.2 httpuv_1.6.15
[45] Wrench_1.20.0 paletteer_1.6.0
[47] magrittr_2.0.3 leaps_3.1
[49] splines_4.3.3 jpeg_0.1-10
[51] doRNG_1.8.6 wk_0.9.1
[53] rootSolve_1.8.2.4 ggbeeswarm_0.7.2
[55] DBI_1.2.2 withr_3.0.0
[57] class_7.3-22 systemfonts_1.0.6
[59] htmlwidgets_1.6.4 fs_1.6.4
[61] SuppDists_1.1-9.7 SingleCellExperiment_1.24.0
[63] ggrepel_0.9.5 labeling_0.4.3
[65] SparseArray_1.2.4 cellranger_1.1.0
[67] flashClust_1.01-2 rncl_0.8.7
[69] see_0.8.4 lmom_3.0
[71] effectsize_0.8.8 zoo_1.8-12
[73] XVector_0.42.0 decontam_1.22.0
[75] secretbase_0.5.0 foreach_1.5.2
[77] fansi_1.0.6 caTools_1.18.2
[79] grid_4.3.3 data.table_1.15.4
[81] ggtree_3.10.1 rhdf5_2.46.1
[83] irlba_2.3.5.1 gridGraphics_0.5-1
[85] DescTools_0.99.54 base64url_1.4
[87] lazyeval_0.2.2 yaml_2.3.8
[89] survival_3.6-4 crayon_1.5.3
[91] RColorBrewer_1.1-3 tidyr_1.3.1
[93] later_1.3.2 codetools_0.2-19
[95] base64enc_0.1-3 GlobalOptions_0.1.2
[97] shape_1.4.6.1 limma_3.58.1
[99] estimability_1.5.1 Rsamtools_2.18.0
[101] ggside_0.3.1 foreign_0.8-86
[103] pkgconfig_2.0.3 ggpubr_0.6.0
[105] xml2_1.3.6 GenomicAlignments_1.38.2
[107] aplot_0.2.2 scatterplot3d_0.3-44
[109] ape_5.8 viridisLite_0.4.2
[111] xtable_1.8-4 interp_1.1-6
[113] car_3.1-2 highr_0.10
[115] hwriter_1.3.2.1 plyr_1.8.9
[117] httr_1.4.7 rbibutils_2.2.16
[119] tools_4.3.3 broom_1.0.5
[121] beeswarm_0.4.0 htmlTable_2.4.2
[123] checkmate_2.3.1 nlme_3.1-164
[125] PCAtest_0.0.1 lme4_1.1-35.3
[127] digest_0.6.36 numDeriv_2016.8-1.1
[129] adephylo_1.1-16 Matrix_1.6-5
[131] farver_2.1.2 reshape2_1.4.4
[133] yulab.utils_0.1.4 viridis_0.6.5
[135] DirichletMultinomial_1.44.0 rpart_4.1.23
[137] glue_1.7.0 cachem_1.0.8
[139] Hmisc_5.1-2 generics_0.1.3
[141] Biostrings_2.70.3 classInt_0.4-10
[143] mvtnorm_1.2-4 spdep_1.3-3
[145] metagenomeSeq_1.43.0 statmod_1.5.0
[147] biomformat_1.30.0 ScaledMatrix_1.10.0
[149] carData_3.0-5 minqa_1.2.6
[151] ANCOMBC_2.4.0 glmnet_4.1-8
[153] utf8_1.2.4 seqinr_4.2-36
[155] gtools_3.9.5 readxl_1.4.3
[157] datawizard_0.10.0 ggsignif_0.6.4
[159] gridExtra_2.3 shiny_1.8.1.1
[161] GenomeInfoDbData_1.2.11 energy_1.7-11
[163] rhdf5filters_1.14.1 parameters_0.21.6
[165] RCurl_1.98-1.14 memoise_2.0.1
[167] rmarkdown_2.26 scales_1.3.0
[169] gld_2.6.6 svglite_2.1.3
[171] phylobase_0.8.12 rstudioapi_0.16.0
[173] cluster_2.1.6 hms_1.1.3
[175] ggcorrplot_0.1.4.1 munsell_0.5.1
[177] colorspace_2.1-0 rlang_1.1.4
[179] s2_1.1.6 DelayedMatrixStats_1.24.0
[181] sparseMatrixStats_1.14.0 circlize_0.4.16
[183] scuttle_1.12.0 mgcv_1.9-1
[185] xfun_0.43 multcompView_0.1-10
[187] coda_0.19-4.1 e1071_1.7-14
[189] TH.data_1.1-2 plotROC_2.3.1
[191] iterators_1.0.14 emmeans_1.10.1
[193] abind_1.4-5 tibble_3.2.1
[195] treeio_1.26.0 gmp_0.7-4
[197] Rhdf5lib_1.24.2 DECIPHER_2.30.0
[199] bitops_1.0-7 Rdpack_2.6
[201] ps_1.7.6 promises_1.3.0
[203] RSQLite_2.3.6 sandwich_3.1-0
[205] DelayedArray_0.28.0 proxy_0.4-27
[207] Rmpfr_0.9-5 compiler_4.3.3
[209] statsExpressions_1.5.4 prettyunits_1.2.0
[211] boot_1.3-30 beachmat_2.18.1
[213] BiocSingular_1.18.0 units_0.8-5
[215] MASS_7.3-60.0.1 kSamples_1.2-10
[217] progress_1.2.3 uuid_1.2-0
[219] BiocParallel_1.36.0 insight_0.19.11
[221] R6_2.5.1 rstatix_0.7.2
[223] fastmap_1.1.1 multcomp_1.4-25
[225] vipor_0.4.7 rsvd_1.0.5
[227] ggstatsplot_0.12.3 nnet_7.3-19
[229] gtable_0.3.5 KernSmooth_2.23-22
[231] latticeExtra_0.6-30 deldir_2.0-4
[233] htmltools_0.5.8.1 RcppParallel_5.1.7
[235] bit64_4.0.5 lifecycle_1.0.4
[237] processx_3.8.4 spData_2.3.0
[239] nloptr_2.0.3 callr_3.7.6
[241] vctrs_0.6.5 zeallot_0.1.0
[243] ggfun_0.1.4 sp_2.1-4
[245] pillar_1.9.0 gplots_3.1.3.1
[247] prismatic_1.1.2 ShortRead_1.60.0
[249] locfit_1.5-9.9 jsonlite_1.8.8
[251] expm_0.999-9 GetoptLong_1.0.5